%%
% 本程序求解静电场问题
%%
close,clc,clear;
% 材料参数设置
epsilon.air = 1;  %---介电常数
%% 
% 激励大小设置(V)
voltage_airout = 0;
voltage_cable1 = 1;
voltage_cable2 = 0;
voltage_cable3 = 0;
%% 
% 读取网格信息
data = load("mesh\T1.mat");  %---导入GmshInf网格文件
coord = data.coord;                     %---坐标信息
coord(:, 2:end) = coord(:, 2:end)*1e-3;  %---坐标单位变换为国际单位制：m
element_s = data.element_S;             %---获取面单元信息
element_l = data.element_L;             %---获取线单元信息

% 常数设定
PosOfcoordx = 2;    %--结点x坐标列
PosOfcoordy = 3;    %--结点y坐标列
PosOfcoordz = 4;    %--结点z坐标列
PosOfcoord = 1;     %---结点编号列
PosOfelement = 2;   %---单元编号列
PosOfLentity = 1;   %---线单元实体编号列
PosOfSentity = 1;   %---面单元实体编号列
PosEleNum1 = 3;     %---单元中结点局部编号1对应的结点全局编号列
PosEleNum2 = 4;     %---单元中结点局部编号2对应的结点全局编号列
PosEleNum3 = 5;     %---单元中结点局部编号3对应的结点全局编号列
PosLEleNum1 = 3;     %---线单元中结点局部编号1对应的结点全局编号列
PosLEleNum2 = 4;     %---线单元中结点局部编号2对应的结点全局编号列

region.air = 1;  %---空气域实体编号
line.airout = 1;  %---空气域外边界实体编号
line.cable1 = 2;  %---电缆边界1
line.cable2 = 3;  %---电缆边界2
line.cable3 = 4;  %---电缆边界3

% 获取网格基本信息
NumOfcoord = coord(end, PosOfcoord)-coord(1, PosOfcoord)+1;  %---结点数目
NumOfelement = element_s(end, PosOfelement) - element_s(1, PosOfelement)+1; %---单元数目
NumOfLentity = element_l(end, PosOfLentity) - element_l(1, PosOfLentity)+1; %---线实体数目
NumOfSentity = element_s(end, PosOfSentity) - element_s(1, PosOfSentity)+1; %---面实体数目
% 区域单元确定
Numid.air = [find(element_s(:, PosOfSentity)==region.air)];  %---空气域单元编号
Numid.airout = [find(element_l(:, PosOfSentity)==line.airout)];  %---电缆边界1单元编号
Numid.lineCable1 = [find(element_l(:, PosOfSentity)==line.cable1)];  %---电缆边界1单元编号
Numid.lineCable2 = [find(element_l(:, PosOfSentity)==line.cable2)];  %---电缆边界1单元编号
Numid.lineCable3 = [find(element_l(:, PosOfSentity)==line.cable3)];  %---电缆边界1单元编号
% patch('vertices', coord(:, 2:3), 'faces', element_s(:, 3:end),'FaceColor',[1 1 1]);
% axis equal;
%%
% 单元形状函数以及单元面积
K = zeros(NumOfcoord, NumOfcoord);  %初始化刚度矩阵
b = zeros(NumOfcoord, 1);  %初始化右端项
be = [coord(element_s(:, PosEleNum2), PosOfcoordy) - coord(element_s(:, PosEleNum3), PosOfcoordy),...
      coord(element_s(:, PosEleNum3), PosOfcoordy) - coord(element_s(:, PosEleNum1), PosOfcoordy),...
      coord(element_s(:, PosEleNum1), PosOfcoordy) - coord(element_s(:, PosEleNum2), PosOfcoordy)];
ce = [coord(element_s(:, PosEleNum3), PosOfcoordx) - coord(element_s(:, PosEleNum2), PosOfcoordx),...
      coord(element_s(:, PosEleNum1), PosOfcoordx) - coord(element_s(:, PosEleNum3), PosOfcoordx),...
      coord(element_s(:, PosEleNum2), PosOfcoordx) - coord(element_s(:, PosEleNum1), PosOfcoordx)];

area = abs(0.5*(be(:, 1).*ce(:, 2) - be(:, 2).*ce(:, 1)));  %--计算每个单元的面积
%%
% 确定刚度矩阵K以及右端项b
for ii = 1:NumOfelement
    if(ii<=Numid.air(end) && ii>=Numid.air(1))
        a_x = epsilon.air;
        a_y = epsilon.air;
        f = 0;
    end
    K(element_s(ii, PosEleNum1:PosEleNum3), element_s(ii, PosEleNum1:PosEleNum3)) = K(element_s(ii, PosEleNum1:PosEleNum3), element_s(ii, PosEleNum1:PosEleNum3))+...
                                                         (a_x*be(ii, :)'*be(ii, :)+a_y*ce(ii, :)'*ce(ii, :)) / (4*area(ii));
    b(element_s(ii, PosEleNum1:PosEleNum3), 1) = b(element_s(ii, PosEleNum1:PosEleNum3), 1)+area(ii)*f / 3;
end
%%
% 添加边界条件
element_airout = unique( element_l(Numid.airout, PosLEleNum1:PosLEleNum2) );
K(element_airout, element_airout) = 1e20*eye(length(element_airout));
b(element_airout) = voltage_airout*1e20;
element_cable1 = unique( element_l(Numid.lineCable1, PosLEleNum1:PosLEleNum2) );
K(element_cable1, element_cable1) = 1e20*eye(length(element_cable1));
b(element_cable1) = voltage_cable1*1e20;
element_cable2 = unique( element_l(Numid.lineCable2, PosLEleNum1:PosLEleNum2) );
K(element_cable2, element_cable2) = 1e20*eye(length(element_cable2));
b(element_cable2) = voltage_cable2*1e20;
element_cable3 = unique( element_l(Numid.lineCable3, PosLEleNum1:PosLEleNum2) );
K(element_cable3, element_cable3) = 1e20*eye(length(element_cable3));
b(element_cable3) = voltage_cable3*1e20;
%%
% 求解
varPhi = K\b;
%%
% 结果云图
patch('vertices', coord(:, PosOfcoordx:PosOfcoordy), 'faces', element_s(:, PosEleNum1:PosEleNum3), 'facevertexCdata', varPhi, 'facecolor', 'interp')
colormap(jet)
colorbar; 
c = colorbar; % 获取颜色条对象
c.Label.String = 'V'; % 设置颜色条的标签
axis equal;
axis image;
xlabel("m");
ylabel("m");
